clear; clc; close all;


%% 1. 参数设置
mu = -0.6;
x0 = 0.8;
m_vals = [-2, -1, 0, 1];  % 对应不同的 q0 = -0.8 + 2*m*pi
k_start = 1.3;
k_end   = 1.8;
dk      = 0.001;
k_vals  = k_start : dk : k_end;

% 迭代相关
transient_steps = 300;   % 舍去的初始过渡步数
record_steps    = 100;   % 每个 k 下记录的迭代点数

% 颜色或标记设置
color_list = {'b','r','m','k'};  % 蓝、红、品红、黑

hold on; grid on;


%% 2. 循环 m_vals，绘制不同初始 q0 的分岔图
for idx = 1:length(m_vals)
    m = m_vals(idx);
    q0 = -0.8 + 2*m*pi;   % memristor 初始状态

    % 用于存放该 m 值下的所有 (k, q) 分岔数据
    k_bif = [];
    q_bif = [];

    for k = k_vals
        % 从 (x0, q0) 出发，先做 transient_steps 步以排除瞬态
        x = x0;
        q = q0;
        for t = 1:transient_steps
            [x, q] = MTM(x, q, mu, k);
        end

        % 再记录 record_steps 步，用于分岔图
        q_temp = zeros(1, record_steps);
        x_temp = zeros(1, record_steps);
        x_temp(1) = x; 
        q_temp(1) = q;
        for t = 2:record_steps
            [x_temp(t), q_temp(t)] = MTM(x_temp(t-1), q_temp(t-1), mu, k);
        end

        % 将 record_steps 个 q 值与同一个 k 组合
        k_bif = [k_bif, k*ones(1, record_steps)];
        q_bif = [q_bif, q_temp];
    end

    % 将本组 (k_bif, q_bif) 以散点或点线绘制
    scatter(k_bif, q_bif, 1, color_list{idx}, 'filled',...
        'DisplayName',['q_0=-0.8 + 2(',num2str(m),')\pi']);


    
    %% 3. 美化图形
    xlabel('k');
    ylabel('q');
    legend('Location','best');
    title('\mu=-0.6, x_0=0.8, q_0=-0.8+2m\pi (m=-2,-1,0,1)');
end

%% ========== 迭代方程：2D-MTM ==========
function [x_next, q_next] = MTM(x, q, mu, k)
    % 2D-MTM (分段映射):
    % if x < 0.5:
    %   x_{n+1} = mu*x + k*cos(q)*x
    % else:
    %   x_{n+1} = mu*(1 - x) + k*cos(q)*x
    % q_{n+1} = q + x
    if x < 0.5
        x_next = mu*x + k*cos(q)*x;
    else
        x_next = mu*(1 - x) + k*cos(q)*x;
    end
    q_next = q + x;
end
